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ABSTRACT 

Unopposed radiative cooling in clusters of galaxies results in excessive mass de- 
position rates onto the central brightest cluster galaxy. However, the cool cores of 
galaxy clusters are continuously heated by thermal conduction and turbulent heat 
diffusion due to minor mergers or the galaxies orbiting the cluster center. These pro- 
cess can either reduce the energy requirements for AGN heating of cool cores, or 
they can prevent overcooling altogether. We perform 3D MHD simulations including 
field-aligned thermal conduction and self-gravitating particles to model this in detail. 
Turbulence is not confined to the wakes of galaxies but is instead volume-filling, due 
to the excitation of large-scale g-modes. We systematically probe the parameter space 
of galaxy masses and numbers to assess when the cooling catastrophe is prevented. 
For a wide range of observationally motivated galaxy parameters, we find that the 
magnetic field is randomized by stirring motions, restoring the conductive heat flow 
to the core. The cooling catastrophe either does not occur or it is sufficiently delayed 
to allow the cluster to experience a major merger that could reset the conditions in the 
intracluster medium. Whilst dissipation of turbulent motions (and hence dynamical 
friction heating) is negligible as a heat source, turbulent heat diffusion is extremely im- 
portant; it predominates in the cluster center. However, thermal conduction becomes 
important at larger radii, and simulations without thermal conduction suffer a cooling 
catastrophe. Conduction is important both as a heat source and to reduce stabiliz- 
ing buoyancy forces, enabling more efficient diffusion. Turbulence enables conduction, 
and conduction enables turbulence. In these simulations, the gas vorticity — which is 
a good indicator of trapped g-modes-increases with time. The vorticity growth is 
approximately mirrored by the growth of the magnetic field, which is amplified by 
turbulence. 

Key words: conduction - cooling flows - galaxies: clusters: general - galaxies: active 
- instabilities - X-rays: galaxies: clusters 



1 INTRODUCTION 

The intracluster medium (ICM) in many galaxy clusters 
has central cooling times shorter then the Hubble time. 
Radiative cooling should lead to large accumulation of cold 
material in their centers; however, there is no observational 
evidence for such gas. This can be understood if some 
source of heating balances cooling in the ICM. The heating 
mechanisms invoked to explain this ove rcooling problem 
involv e AGN "radio mod e" heat ing (e.g.. iBinnev fc Taborl 



iRuszkowski et all l|2004; 



19951); Ichurazov et all |2002l ); iFabian et al l J2005t): 



Scannapieco fc Briiggcn 



1 20081 )). preheating by AGN ^McCarthy etall |200S|), cos - 
mic rays from AGN l|Guo fc O h 2008; Sharma et al.ll2009al ). 



supernovae, turbulent mixing dKim fc Naravanl l2003al : 
IVoigt fc Fabianl 12004 iDennis fc Chandranl 120051), thermal 
condu ction |Zakamska fc Naravanl 120031 : iKim fc Naravanl 
2003b [), a combination o f ther mal conduction and AGN 
( Rus zkowski fc Begelman] 20021) and d ynam i cal fr iction 
fel-Zant et alj 120041: IKim et al.l 120051: iKiml 120071); see 



IConrov fc Ostrikerl (|200SD and lMcNamara fc Nulsenl (|2007l ) 
and references therein for reviews of the above mecha- 
nisms). 

Conduction alone is unlikely to offer the complete 
solution to the overcooling problem for the full range 
of cluster masses, as its strong temperature dependence 
implies that it is less effective in lower mass clusters. 
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Furthermore, thermal conduction is well known to be 
an unstable heating mechanism, either failing to avert a 
cooling catastr ophe, or leading to an is othermal temper- 
ature profile (jBregman fc Davidl Il988l ; iGuo fc Ohl 120081; 
IConrov fc Ostriker 20081 ). Nevertheless, thermal conduction 
may entirely suppress cooling in non cool-core (NCC) 
clusters and reduce the constraints on the requi red energy 
inject ion by AGN in cool-core (CC) clusters (|Guo et al.l 
I2008T ): indeed, it is differe nt to stabilize massive cl usters 
with AGN feedback alone l|Conrov fc Ostriker! liooH ) . and 
a second heat source (such as conduction) is generally 
required. Besides offsetting radiative losses and stemming 
a cooling catastrophe, conduction can have important 
implications for esta blishing the obse rved bimodality in 
cluster core entropy (|Guo et alj [2008) , and th e star for- 
matio n threshold in brightest cluster galaxies (|Voit et al.l 
2008). Indeed, a sudden increase in conduction (due to 
say, turbulence from an AGN outburst or a me rger) could 
mediate a cool core to non cool-core transition l|Guo fc Ohl 
l2009h . 

Thermal conduction can be strongly suppressed by 
magnetic fields that are known to be present in the fCM 
(e.g., Enfilin et al. (2003), Vogt et al. 2003). However, inter- 
est in thermal conduction as a potential hea ting mechanism 
was revived bv lNaravan fc Medvedevl (j2001) who suggested 
that even in the presence of tangled B-fields, the level of 
conduction can be an appreciable fraction of the Spitzer- 
Braginskii value. Computing the exact magnitude and dis- 
tribution of the effective conductivity of the 1CM is fur- 
ther complicated by buoyancy instabilities which re-orient 
the magnetic field. When temperature increases in the di- 
rection of gravity, as in th e cluster outs k irts, the magne- 
tother mal instability (MTI; iBalbusl (l200d ); |Parrish fc Stone] 
(2005)) tends to make the B-fields radial, thereby increasing 
the effective conduction. On the other hand, in cool cores 
where temperature decreases in the dir ection of gravity , 
the heat flux buoyancy instability (HB I; Quataerj (2008); 
IParrish et all (120091 ); iBogdanovic et all ||2009l )) tends to re- 
orient the fields in the direction perpendicular to that of 
gravity, effectively shutting down thermal conduction. 

However, these instabilities do not operate in a static 
atmosphere. Chandra and XMM observations show that the 
cluster gas is rarely in perfect hydrostatic equilibrium. Slosh- 
ing motions due to minor mergers, AGN, or galaxy motions 
can continuously and significantly perturb the g has 
been re peatedly seen in many disparate numerical simu- 
lations tevrardlll990l; iNorman fc Brvanl Il999l; iNagai et all 
2003l;IVazza et al. 2009, 2010l;lAscasibar fc Markcvitch 2006; 



ZuHone et al.l [2010). Current observational evidence for 
turbulence ranges fro m the analysis of pressure maps 
l|Schuecker et alj 12004), its effect on resonant-line scatter- 
ing (Ch urazov et al. I |2004|). and Faraday rotation maps 
I Vogt fc EnfllinlbOOa ; Enfflin fc Vogtll2006l). as well a s con- 
straints on turbulent line widths (|Sanders et ai1 l2009). Low 
levels of turbulence in the ICM can randomize the field con- 
figuratio n set up by the HBI and restore the heat flow to 
the core (|Ruszkowski fc Ohll2010l ; |Parrish et alj|2010l ). Both 
of these works modeled the ICM turbulence via a simple 
driving mechanism to determine the level of turbulence re- 
quired to effectively restore thermal conduction. This ap- 
proach did not allow us to link the level of turbulence to 
the physical properties of the cluster (such as mechanical 



luminosity of the central AGN or the properties of clus- 
ter galaxies). Also, the driving mechanism led, by construc- 
tion, to volume-filling turbulence which was very effective 
in randomizing the magnetic field. While the low ampli- 
tude of the required subsonic turbulence is eminently fea- 
sible (vt ~ 150km s _1 ~ 0.1c s ), the realism of volume-filling 
turbul ence is less clear. For inst ance, both analytic calcula- 
tions ( Subramanian ct aL[|2006j) and numerical simulations 
l|lapichino fc Niemeverll2008l ) predict that turbulence due to 
galaxy wakes should not be volume-filling (/v Ss 0.2 — 0.3), 
as turbulence is largely confined to 'streaks' behind orbiting 
galaxies. 

In this paper, we extend our previous work and per- 
form three-dimensional MHD simulations of the effect of 
turbulence driven by galaxy motions on the properties of the 
anisotropic thermal conduction. We show how the trapping 
of gravity modes excited by the orbiting galaxies can lead to 
volume-filling turbulence of the right magnitude to restore 
conductive heat flow. We demonstrate how these subsonic 
motions generate vorticity and lead to the growth of mag- 
netic field via kinematic dynamo action. We also show that 
turbulent heat diffusion is an important part of the energy 
budget. 

The paper is organized as follows. In [J2]we review basic 
theoretical expectations for the interaction between turbu- 
lence and magnetic fields. In |3]we describe the numerical 
methods and the setup of the inital conditions. In |4] we 
describe our results, including the level and volume-filling 
nature of turbulence, evolution of the gas temperature, 
generation of vorticity and magnetic fields, and nature of 
heating mechanisms. Conclusions are summarized in fj5] 



2 THEORETICAL EXPECTATIONS: 
TURBULENCE AND TRAPPING OF 
G-MODES 

It is useful to begin by reviewing some basic theoretical ex- 
pectations for the behavior of turbulence excited in galactic 
wakes in clusters. In principle, orbiting galaxies can excite 
galactic wakes by two means: hydrodynamically (as the ICM 
collides with the ISM of the galaxy), and gravitationally 
(similar to dynamical friction for collisionless particles). In 
practice, we shall conservatively assume that ram pressure 
stripping is efficient in removing the ISM of galaxies and 
thus that galaxies only exert gravitational influence. 

Volume-filling Turbulence The volume filling factor 
of galaxy wakes is small ( for a simple analytic estimate, see 
ISubramanian et al.1 i|2006l )). This might seem to imply that 
the impact of turbulence excited by galactic wakes is con- 
fined to a small fraction of the cluster. However, orbiting 
galaxies can also resonantly excite g-modes, which from a 
formal WKBJ analysis have the dispersion relation: 



2 fei 



(1) 



k 2 = k\ + k 2 . The Brunt- Vaisala, frequency for buoyant 
oscillations is: 



hydro \ 2 
J BV ) : 



'> a H 



3dlnS dlnT 



5 dlnr dlnr 



(2) 
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where S is the fluid entropy S = feT/n 2 ' 3 , and oj^y™ ap- 
plies if thermal conduction is negligible, while lub 1 ™ applies 
if the thermal conduction time is sufficiently short that a 
displaced blob's temperature i s determined by con ductive 
rather than adiabatic cooling |Sharma et al.l l2009bl ). Note 
that WBv r °! a; BV ° depend on the entropy and temperature 
gradient respectively; typically, uj^™ is about a factor of 2 
smaller than cj B y lro . 

The above dispersion relation immediately implies that 
to obtain modes with real k r , the driving frequency to < 
wbv; otherwise the modes have imaginary radial wavenum- 
ber and are evanescent. Physically, one can always achieve 
a low frequency response by making the mode progres- 
sively more tangential, but it is impossible to drive the sys- 
tem at frequencies higher than the maximum response fre- 
quency of cjbv, corresponding to completely vertical oscil- 
lations. This thus implies that waves driven at frequencies 
cj < cjbv can be resonantly excited, and must propagate 
inward toward the cluster cente r (as c an be seen from their 
group velocity; [Balbus fe Sokeil (|l990h ). where they will be 
trapped, reflected and focused insid e the resonance radiu s 
where ljbv = w. A linear analysis bv lBalbus fc Sokerl (| 19901 1 
showed that most of the power in (/-modes is in the longest 
wavelength^]. Note that both ui (which depends on the or- 
bital frequencies of galaxies) and ljbv are sensitive to the 
gravitational potential, which is instrumental in determin- 
ing if g-modes will be excited. 

Isotropic Turbulence Turbulence in the fluid has to 
compete with buoyancy forces arising from stable strati- 
fication. One can show that the ratio of tangential and 
radial velocities is given by (e.g., see discussion in §2 of 
iRuszkowski fc Ohl (|2010h ): 



— ~ ( CUl 



V^BV 



Fr" 



(3) 



where cjl = v/L is the eddy turnover frequency at a given 
scale, and Fr is the Froude number, which compares inertial 
and gravitational forces (Ri ~ 1 /Fr 2 is the Richardson num- 
ber). If uj <C wbv, then turbulence is fundamentally 2D, and 
for instance it is difficult to rearrange magnetic fields in the 
radial direction. However, the level of turbulence required 
to overcome stable stratification is weak; for typical cluster 
condit ions the critical turbulent velocity is l)Sharma et al.l 
l2009bl ): 



a « 135 kms <?_g 



1/2 1/2 



dlnT/dlnr 
015 



1/2 



Ri \" 1/2 

— I (4) 
0.25 J v ' 



where g_g is the gravitational acceleration in units of 
1CP 8 cm 2 s -1 , no is a characteristic scale height in units of 10 
kpc, and Ri c is the critical Richardson number; Ri c ~ 1/4 
is typical for hydrodynamic flow. 

At first blush, the requirement for Fr > 1 might seem 
to be at odds with the requirement that lv < ljbv for g 
modes to be excited. However, note that for homogeneous 
Kolmogorov turbulence, cjl oc L~ 2//;! ; it is therefore con- 
ceivable that low frequency g-modes can be excited on large 



1 Although a WKBJ analysis form ally breaks down in this 
regime, a subsequent numerical study jLufkin et al I 1995ft showed 
that many of the linear theory results are still valid. 



scales, while high-frequency small-scale modes can overcome 
stabilizing buoyancy forces. Since our background state is 
not homogeneous, we have to resort to 3D simulations to 
verify if this expectation is indeed satisfied. This is a major 
goal of this paper. 

Vorticity and B-field growth G-modes excite vor- 
ticity. An easy way to see this is to examine the vorticity 
evolution form of the mo mentum equation fo r g- waves (i.e., 
assuming 8P/P < Sp/p) (|Lufkin et alj|l995h : 



(5) 



where fi is the vorticity, and to note that k is in general non- 
radial, so that k x g is non-zero (indeed, we see in Fig. [2]that 
since lj/lobv rises toward the center, that g-modes become 
progressively more tangentially biased there). This implies 
that vorticity is a good tracer of g-modes, a fact which we 
shall exploit. It also means that g-modes could conceivably 
drive an efficient dynamo. There is a well-known analogy 
between the vorticity equation: 

— + V x (fi x u) = -V x (z/V x n) (6) 

where v is the viscosity, and the relation for the magnetic 
field in the flux- freezing limit, 

<9B 



dt 



+ V x (B x u) = -V x (77V x B) 



(7) 



where r] is the electricity resistivity. This, together with the 
fact that the divergence of B and fl both vanish, leads to 
the expectation that their growth might be relatecQ There 
have been a number of studies pointing out that turbu- 
lent motions could giv e rise to magnetic fields in c luster s 
(e.g.. [ Ruzmaikin et al.l dl989T ) ; [ Subramanian et al.l |2006); 
iRvu et all |200Sl '); ICho et all (|2009l ')'). This subject is rich 
and beyond the scope of this paper; we shall merely com- 
pare the growth of vorticity and magnetic fields in our simu- 
lations, to see how well they track one another. A reasonable 
expectation is that the magnetic fields achieve equipar tition 
with turbulence (e.g., ISchekochihin fc Cowlevl ((2007) , and 
references therein): 



^turb 



100 kms- 



(8) 



where «t ur b is the rms turbulent velocity on large scales. 
The above estimate is consistent with observed ~ pG fields 
(|Carilli fc Tavlorl I2002T I , though there are considerable un- 
certainties The fact that trapping of g-modes can give rise 



2 Note, however, that this analogy is imperfect, since fl = V x u, 
which leads to a nonlinear coupling in the equations, whereas no 
such relation exists between B and u. 

3 In general, estimates based on rotation measure (RM) lead to 
stronger magnetic fields, while those based on synchrotron and 
inverse Compton (IC) analysis give weaker fields. However, RM 
methods may overestimate fields if single-scale magnetic field cor- 
relation length is used llNewman et alj|2002ft or when the small- 
scale fluctuations in den sity and magnet ic field are correlated 
in a turbulent medium llBeck et all [20031. Moreover, these es- 
timates depend on whether radio sources used to probe the field 
strength are embeded in the ICM, with smaller values inferred 
when background sour ces rather than embedded ones are used 
dCarilli fc Tavlodl2002ri . 
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to volume-filling turbulence would then be instrumental in 
allowing volume-filling magnetic fields. 

Mag netic tens ion Magnetic tension can inhibit the 
HBI (Quatacrt 2008). For perturbation scales comparable 
to the radius r (i.e., A = 2r) we obtain a critical value: 



lO^G, ( 



10 "cms 

1/2 



1/2 



30 kpc 



/dlnT/dlnr 
V 03 



0.02 cm" 3 

1/2 



1/2 



(9) 



for suppression, where the fiducial values are measured from 
our simulated cluster at a radius of r = 30kpc. Due to the 
similarity between the field values in equation @ and ©, it 
has been suggested that magnetic fields amplified by turbu- 
le nce can prev e nt the onset of the HBI (e.g., see discussion 
in iKunz et al.l j2010f ). Note that their version of equation 
([9]) yields somewhat lower B-fields than ours, for identical 
parameters. In any case, equation (£9|) is only an approxi- 
mation as the derivation assumes the WKB approximation, 
while the non-linear saturation of the HBI occurs on global 
scales). In this paper, we will deliberately ignore this pos- 
sibility. Observationally, the strength of the magnetic field 
in the ICM is ~ /iG and has a large scatter of about an 
order of magnitude wi thin the ICM and between clusters 
(jCarilli fc Tavlorl I2002T ); moreover, there are considerable 
observational uncertainties in these values, as mentioned 
above. Numerical simulations show that the HBI still devel- 
ops for ~ /iG fields (I. Parrish, priv. comm.), although it can 
be delayed for increased field strengths. Given the large un- 
certainty in whether observed field strengths are capable of 
stabili z ing the HBI, past s tudies of HBI (e.g.. iParrish et aU 
j2009h;lBogdanovic et ail (|2009l ); iRuszkowski fc Obi (|201(t ): 
IParrish et alj l|2010] n focused on the regime where the mag- 
netic tension is unimportant. We also adopt the same ap- 
proach here, and study if volume-filling turbulence alone 
can stabilize the HBI. More specifically, we consider plasma 
/) > 1 and note that, as long as the field is not dynami- 
cally important, its exact value does not play a role. In this 
case, the magnetic field strength scales out of the problem 
and only serves as a medium to redirect the heat flow via 
anisotropic thermal conduction. We can therefore study the 
effects of turbulence alone without the possibly confounding 
effects of magnetic tension. 

Turbulent heating and heat diffusion Turbulence 
impacts the thermodynamics of the fluid through its effect 
on thermal conduction, both randomizing and amplifying 
the magnetic field. Both of these suppress the HBI, and al- 
low thermal conduction at ~ 1/3 the Spitzer value. How- 
ever, turbulence can also directly affect the thermal state of 
the plasma through dissipation of turbulent motions (direct 
heating), or allow i ng hea t transport via turbulent diff usion 
|Kim fc Naravanl (|2003bl ): iDennis fc Chandranl (|2005T ). and 
references therein). The heating rate from dissipation of ki- 
netic and magnetic energy is: 



Tdis; 



I 



(10) 



where Cdi ss is a dimensionless constant of order unity and I, 
the dominant velocity length-scale, is unknown but almost 
certainly a fu nction of radius; a reasonab le ansatz might be 
I ~ ctr + lo jDennis fc Chan dran 2005J), where a is some 
adjustable constant of order unity, and lo is some minimal 



lengthscale. On the other hand, the heating rate from tur- 
bulent heat diffusion is: 

r diff = V • (KturbpTVs) (11) 

where s = Cvhi(p/p 7 ) is the specific entropy, and the tur- 
bulent diffusivity is: 



Kturb ~ ul min (1, ( — ) ) 

V^bv, 



(12) 



where the second factor of (lo/ljbv) 2 takes into account 
the damping of radial hea t transport by buoyancy forces 
(Den nis fc Chandranl [2005h . The fact that K tur b ~ ul is 
of order the hydrodynamic value even for a magnetize d 
plasma was found in MHD simulations bv lCho et all (2003). 
Nonetheless, equation (|ll|l should be understood to be 
only approximate, since it assumes that fluid elements 
are transported adiabatically, which need not be the case 
when anisotropic conduction is operating. In reality, both 
the thermal conduction diffusion coefficient respitzcr = 
Velmtp ~ 10 3 °cm 2 s _1 (n/10~ 2 cm- 3 )- 1 (T/2keV) 5/2 
and the turbulent heat diffusion coefficient Kturb ~ 
10 30 cm 2 s _1 (it/200kms _1 )(//20kpc) can be comparable, 
and either could dominate in a specific situation. 

Thermal conduction may indirectly assist with tur- 
bulent heat diffusion, as it reduces the impact of buoy- 
ancy forces (and thu s reduces ubv)- Indeed, simulations by 
ISharma et all (2009) show that metal mixing in a stratified 
plasma is much more efficient once conduction is at play, 
allowing much broader metallicity profiles, for this very rea- 
son. Naively, if we think of gas entropy as a scalar to be 
advected by turbulent motions, similar conclusions should 
hold, although of course the interaction between heat trans- 
port and dynamics requires detailed simulations. We shall 
investigate the relative role of all these heating processes in 
our simulations. 



3 METHODS 

3.1 Initial conditions for the gas 

The details of the numerical setup are described in 
Ruszkowski & Oh (2010; hereafter RO10). Here we summa- 
rize key differences. The cluster parameters used here are 
similar to those corresponding to cool-core cluster A2199. 
In addition to the NFW potential of the cluster halo, we 
also include the contribution from the central brightest 
cluster galaxy (BCG), which was not included in RO10. 
The gravitational potential is described by the sum of the 
term due to an NFW profile with a softened core 



-2GM 



^ 1 + r/r c ln(l + r/r c ) 



1 + r/r s 



-2GM Q 



r s (r s — 2r c ) ln(l + r/r s ) 



rAr s — r c 



r/r c 



(13) 



where r c is the smoothing core radius (r c = 20 kpc), 
r s = 390 kpc is the usual NFW scale radius, and the BCG 
contribution which has a King profile: 



-9cr 2 c 



In (a: + y/l + x 2 ) 



(14) 
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where x = r/r^ cg , r^ cg = 3 kpc is the core radius for the 
BCG and (Jbcg = 200 km s _1 is its line-of-sight velocity 
dispersion. The parameter M = 3.8 x 1O 14 M0 in equation 
(|14[) determines the cluster mass and is of the order of the 
total cluster mass, M200 = 6.6 x 10 14 M§. We then solve the 
equation of hydrostatic equilib rium assuming the en tropy 
distribution as parametrized bv lCavagnolo et al.l (|2009l ); see 
equations (16) & (17) of RO10. Note that we do not include 
the gravitational contribution from other galaxies ( £I3,2[I in 
our initial conditions, so the system is not initially in full 
hydrostatic equilibrium. However, after an initial transient, 
it rapidly relaxes to a new equilibrium configuration. 

The addition of the BCG has two effects. Firstly, due to 
the increased gravitational acceleration, it results in higher 
gas densities compared to the models we considered in 
Ruszkowski & Oh (2010). This allows for a more conser- 
vative analysis of the effect of cooling. In fact, the central 
density here is a factor of ~ 3.5 times higher, which, com- 
bined with a slightly lower assumed central temperature, 
results in a central cooling time which is nearly 5 times 
shorter. The higher adopted central density in this paper i s 
in line with that observed in A2199 ll Johnstone et al.1 l2002). 
Given this more stringent setup, some of the stable models 
in Ruszkowski & Oh (2010) would actually undergo a cool- 
ing catastrophe. Secondly, the change in the gravitational 
potential has consequences for the Brunt- Vaisala frequency 
and the trapping of g-modes, as we discuss below. 

The initial distribution of density and temperature is 
shown in Figure [T] The frequency of circular orbits cj or b 
and the Brunt- Vaisala frequencies uj^^ t °, for a hy- 

drodynamic and conducting fluid with this initial density 
and temperature profile are shown in Figure [2] For a mode 
with a given value of cj, g-modes can be resonantly excited 
if ui < ojbv- This therefore defines an outer trapping radius 
for such a mode. Note that both uj OT h and ujbv are both 
strong functions of the gravitational potential. We have di- 
rectly verified the resonance condition by running simula- 
tions both with and without the central cD galaxy; in the lat- 
ter case, orbiting galaxies fail to excite volume-filling turbu- 
lence, which is to be expected since ujbv falls inward in this 
case, and the reson a nce conditio n is never satisfied (see also 
iLufkin et ail l| 19951 ) i lKiml (120071 )1. Note that fine-turning of 
the resonan ce condition is not nece ssary: the resonance is not 
very sharp (|Balbus fc Sokerlll990h . and in practice galaxies 
with non-circular orbits excite modes with a variety of har- 
monics, some of which can potentially fall below cjbv- 

The magnetic field setup was identical to that in 
Ruszkowski & Oh (2010): we generate statistically isotropic 
random-phase complex fields in Fourier space, with 3D 
Fourier amplitudes given by: 



0.1000 F ' ' 1 ^H10 



E 




0.0001 I i i i . 1 

1 10 100 

radius [kpc] 

Figure 1. Initial electron number density (solid line) and tem- 
perature (dashed line) in the ICM of our simulated cluster. 




0.1 1.0 10.0 100.0 1000.0 

radius [kpc] 



Figure 2. The frequency of circular orbits cj or b (solid line), the 
Brunt-Vaisala frequency lJ^ to (dashed line) for a hydrodynamic 
fluid and cj^y D (dotted line) for a magnetized conducting fluid 
(all in [Hz]). The frequencies correspond to the initial density and 
temperature profile shown in Figure 1. 




as appropriate for Kolmogorov turbulence, where k = 
2tt/\ and A ~ 43/i -1 kpc is a smoothing wavelength. We 
then apply a divergence cleaning operator in k— space, and 
then inverse Fourier transform the field back to real space. 



3.2 Initial conditions for the galaxies 

The simulations must be initialized with a galaxy popula- 
tion, which has the appropriate spatial distribution, masses 
and velocities. Rather than relying upon cosmological simu- 
lations, we use an empirically grounded approach, which also 
has the advantage of speed and flexibility. How are the galax- 
ies spatially distributed? From a samp le of K-band sel ected 
galaxies within 93 clusters and groups. [Lin et all l|2004h find 
that the galaxy number density profile in clusters is well 
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described by the NFW profile l|Navarro et alJll997h with a 
concentration parameter c ~ 3, with no evidence for cluster 
mass dependence of the concentration. The theoretical jus- 
tification for galaxies tracing the NFW profile is somewhat 
equivocal. If one attempts to use cosmological simulations 
to set up initial conditions, the radial distribution of sub- 
halos in simulations is well-known to b e less concentrated 
than the dark-matter, or 'anti-biased' l|Nagai fc Kravtsovl 
(2005), and references therein). This is due to tidal stripping 
of sub-halos in the central regions. On the other hand, sim- 
ulations that include galaxy formation allow subhalos to be 
selected by stellar mass. This generally shows closer agree- 
ment with observed profiles, as the stellar mass (which is 
tightly bound) remains con served while the dark m atter is 
stripped from outer regions (|Nagai fc Kravtsovll2005l ). Other 
investigators find that the fraction of such stellar-dominated 
halos is small, but caution that numerical resol ution effects 
may preclude robust conclusions at this point (|Dolag et al.l 
2009). Overall, we therefore simply employ the observational 
result that galaxies trace the NFW profile. 

As for the galaxy masses, instead of using a 
Schechter function, we sim p ly assume ( as did, for instance, 
ISubramanian et al.l (120061 ); iKiml i|2007l )) that all galaxies 
have the same mass. This is for two reasons. Firstly, this al- 
lows us to rapidly explore the effect of varying galaxy masses 
(due, for instance, to different efficiencies of tidal stripping). 
The assumption of a characteristic mass is reasonable: since 
dynamical friction scales as M 2 al , most turbulent motions 
are induced by galaxies of mass ~ M*, where most of the 
mass resides, rather than the more abundant lower mass 
galaxies. Indeed, we shall find that the induced gas motions 
are mostly sensitive to the mass of galaxies, and less sensi- 
tive to their number ( H4.ip . Previous hydrodynamic simu- 
lations found unchanged results with galaxies drawn from 
a Schechter distribution , if the characteristic break mass 
M* ~ M ga i (|Kimll2007T l. Secondly, it allows us to directly 
calibrate against lensing estimates for subhalo mass frac- 
tion. Unlike K-band surveys, lensing is directly sensitive to 
total mass, but is generally only sensit ive to subhalos with 
M £ 10 11 M Q . iNataraian et al.l ([2009) find from the mas- 
sive lensing cluster CI 0024+16 that ~ 30% of the cluster 
mass can be attributed to substructure with M > 10 11 Mq, 
with typical masses ~ 10 12 Mq (with a weak radial trend 
such that galaxies in the outer regions are more massive; 
see their Fig 6). Their results, including the mass function 
as a function of radius, is broadly consistent with the re- 
sults of the Millenium simulation run, except that the typ- 
ical masses of galaxies is lower in simulations by a factor 
of ~ 3. This is subject to the uncertainties of extra bind- 
ing due to a compact stellar halo mentioned above; note 
that masses of ~ 10 12 Mq is also consist ent with other ob- 
servations from lensin g (IShin et al.ll200sl ) and galaxy wakes 
l|Sakelliou et al.ll2005l ). Below we explore a grid of models 
with varying galaxy mass, but never allowing the total sub- 
structure mass fraction to rise above ~ 25%. For simplicity 
in the code, the galaxies are modeled as point masses. Since 
we are primarily concerned with the excitation of g-modes 
on scales much larger than ~kpc galactic scales, we do not 
expect this simplification to significantly impact our results. 

Given these assumptions, the most rigorous way to ini- 
tialize galaxy velocities is to directly construct the distri- 
bution function from the density profile, using Eddington's 



formula |Binnev fc Tremainell20o3 ; lKazantzidis et al.ll2004h . 
However, velocity anisotropy is only easily incorporated in 
such models if it has certain parametric forms, as for in- 
stance in Osipkov-Merritt models. Instead, we construct a 
self-consistent velocity model via the local Maxwellian ap- 
proximation: approximating the velocity dispersion tensor 
by a multi-variate Gaussian at each point, wi th disper- 
sions given by the solution of the Jeans equation (|Hernquistl 
Il993h . This has the virtue of simplicity and flexibility. Note 
that such models may not be in strict equilibriu m, and 
can demonstrate evolution |Kazantzidis et all 120041 ). How- 
ever, ISpringel et ail (|2005D find the actual amount of relax- 
ation to be small; furthermore, iFaltenbacher et "all l|2005h . 
who directly simulate the motion of galaxies in clusters, 
find their velocity distribution is indeed closely Maxwellian, 
with good agreement between simulation results and equilib- 
rium Jeans equation solutions. We therefore solve the Jeans 
equation assuming no rotational support or bulk streaming 
(v r = V4, = V0 = 0) : 



(n ga i<T 2 ) + 2/3 v 



dc/> 

r dr 



Mgai dr 

where f} v is the velocity anisotropy parameteijf]: 



(16) 



Pv{r) = l 



2<j2(r)' 



(17) 



?i ga i is the galaxy number density and <f) is the combined 
cluster + cD galaxy gravitational potential. Note that we 
have not built self-consistent models and ignore the contri- 
bution of galaxies to the gravitational potential; for a large 
sub-halo mass fraction, the system is not in full equilibrium. 
In practice, this is a small effect, and the galaxy distribution 
does not evolve significantly over the course of our simula- 
tion. 

What are appropriate assumptions for /3„(r)? It may be 
estimated from observations via Jeans equation modelling, 
given knowledge of galaxies positions, the cluster potential, 
and line of sight velocities. A detailed study of 10 clusters 
using a spectroscopic sample of galaxies from SDSS and 
2dF found galax y orbits to be isotr opic within the errors 
for most clusters l|Hwang fc Ledl2008l ). An earlier paper, us- 
ing ENACS data, found that the brightest ellipticals do not 
yield an equilibrium solution, while other ellipticals, SOs 
and early spirals have iso tropic orbits, and late spi rals prefer 
radial to isotropic orbits (|Biviano fc Katgertll20M ). Overall, 
we assume isotropy /3„(r) = 0, and regard this as our default 
model. In passing, we note that one could easily incorporate 
the effects of the velocity anisotropy by, for ex ample, con- 
siderin g fits to to measurements in simulations |Hoeft et al.l 
(2004)). Given that the evidence for orbital anisotropy in 
observations is marginal to date, we defer the study of the ef- 
fect of such orbital distributions to future work. We also note 
that preferentially radial orbits would enhance the restora- 



4 Prom cosmological simulations. iBensonl {2005) finds that radial 
and tangential velocities can be correlated (at least at the time 
of merger), a detail we ignore. 
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tion of conduction even further and strengthen our conclu- 
sions. 

We solve equation (|17[) as an initial value problem, 
where ol{r 200) ~ GM2oo/3r2oo- Having solved for a r (r) and 
Ot(r), we randomly sample from the multi-variate Gaussian 
at each position to create a realization of the velocity field. 
The above procedure allows us to initialize the simulations 
with a realization of galaxies with both masses and six- 
dimensional phase space coordinates (position and velocity). 



3.3 Simulation 

The simulations were performed using the FLASH code 
(version 3.2). FLASH is a modular, parallel adaptive mesh 
refinement magnetohydrodynamic code. Magnetic field 
evolution was solved by means of a directi o nally unsplit 
staggered mesh algorithm (USM; iLee et all l |2009h ). The 
USM module is based on a finite-volume, high-order Go- 
dunov scheme combined with constrained transport method 
(CT). This approach guarantees divergence-free magnetic 
field distribution. We implemented th e anisotropic conduc- 
tion u nit following the approach of ISharma fc Hammettl 
(2007). More specifically, we applied monotonized cen- 
tral (MC) limiter to the conductive fluxes. This method 
ensures that anisotropic conduction does not lead to 
negative temperatures in the presence of steep temperature 
gradients. The three-dimensional computational domain 
was approximately 1 Mpc on each side, enclosing a large 
fraction of the cluster. The central regions of the cluster 
had an enhanced refinement level. The maximum spatial 
resolution for 6 levels of refinement was ~ 2.7h~ kpc. 
The simulations were performed on a 384-processor cluster 
located at the Michigan Academic Computing Center at the 
University of Michigan in Ann Arbor and on the Columbia 
supercomputer at NASA Ames. 



4 RESULTS 

We performed a total of 16 runs including radiative cooling, 
anisotropic thermal conduction and self-gravitating parti- 
cles to emulate the gas "stirring" by galaxies. We consid- 
ered a uniform grid of parameters: 50, 100, 150, and 200 
galaxies characterized by masses of (0.3, 0.6, 0.9, 1.2) xlO 12 
M©. With our cluster mass of 6.6 x 10 14 Mq, these parame- 
ters corresponds to a mass fraction in galaxies ranging from 
/g a i = 2.2% to a maximum of / ga i = 27%. For instance, 
for galaxies with mass 6 x 10 11 Mq, our grid corresponds 
to /gai = (4.3,8.3, 12, 15)%. We also performed two control 
runs: one without the galaxies (and hence without stirring) 
to isolate the effect of heat buoyancy instability, and one 
without conduction to isolate the effect of dynamical fric- 
tion heating by galaxies. 

4.1 Gas velocities and volume filling of turbulence 

Figure [3] (left panel) shows the evolution of the velocity 
dispersion measured within 100 kpc from the cluster cen- 
ter. Thin blue (red) lines are for 100 (200) galaxies respec- 
tively and for equally spaced masses ranging from 3 x 10 11 to 



1.2 x 10 . The mass increases gradually from the lightest to 
the darkest color. The black dashed line is for the pure HBI 
case. The HBI case and lighter-colored curves are evolved 
for shorter times. These runs suffer from overcooling and the 
central temperatures reaches the low temperature threshold 
at which point the simulation is stopped. The right panel 
in Figure [3] shows the median velocity within 100 kpc; the 
color coding corresponds to that in the left panel. It is clear 
from these figures that there is a clear trend for the velocity 
dispersion or the median velocity to increase with the typi- 
cal galaxy mass. A similar, albeit weaker, trend is seen for 
the galaxy number. This is consistent with the findings of 

Kim (2007) who found in pure hydrodynamic simulations, 

1/2 

that the gas velocity dispersion a scales as a oc N ^ A/ ga i, 
where iVgai and M ga i are the number and mass of galaxies, 
respectively. Note that a scaling Ek oc a 2 oc Ng^M 2 ^ is con- 
sistent with dynamical friction in the linear regime, since 
Ek oc M^i for dynamical friction. 

As iVgai and M ga i increase, the cooling catastrophe is 
delayed, and is completely staved off at the upper enve- 
lope of these parameters. In this re spect our M HD simula- 
tions differ markedly from those of iKiml |2007l ). who found 
that a cooling catastrophe was inevitable in purely hydrody- 
namic simulations, for all portions of parameter space. We 
explore these differences further in ij4.5l Note that in our 
case the velocity dispersion seems to increase more slowly 
than a oc N 1 ^ 2 . Besides the inclusion of MHD in our sim- 
ulations, differing results could be due to a variety of fac- 
tors, including the different assumed distribution of galaxies. 
Note that the stated number of galaxies are distributed over 
the entire cluster; the number of galaxies in the inner re- 
gions which actually result in trapped g-modes is actually 
considerably smaller, and subject to Poisson fluctuations. 
Also, the introduction of more galaxies and/or increasing 
their mass does not cause the velocity dispersion to increase 
without limit; inst ead, the gro wth in velocity dispersion ap- 
pears to saturate. IKiml ((2007) also observed this in his hy- 
drodynamic simulations, and attributed it to loss of resonant 
excitation once density fluctuations become large and the 
background is nonlinear. We see the same saturation on the 
same ~ 10 8 yr timescale, in simulations with driven volume- 
filling turbulence, where resonant excitation of modes is not 
an issue (see §3.4 of RO10). The asymptotic velocities of 
~ 100 — 200kms _1 , while generally insufficient for turbu- 
lent heating to be important, is enough to restore thermal 
conduction and enable turbulent heat diffusion. 

The comparison between the velocity dispersion and 
median velocity reveals that that both of these quantities 
are comparable, as might be expected for volume filling tur- 
bulence. 



4.2 Bias in magnetic field orientation 

The evolution of the anisotropy /3 in the orientation of 
magnetic fields is shown in Fig. [4] The definition of this 
parameter is similar to that for galaxy velocity /3 V de- 
fined in Eq. (16). The only difference is that the veloc- 
ity dispersions are replaced by magnetic field dispersions. 
Thus, ft = corresponds to isotropic magnetic fields, whilst 
ft — > (—00, 1) corresponds to progressively more tangential 
(radial) fields respectively. Thin blue lines are for 100 galax- 
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Figure 3. The evolution of the velocity dispersion (left panel) and median velocity within the central 100 kpc (all in [km s^ 1 ]). Blue 
(red) lines are for 100 (200) galaxies respectively, for equally spaced masses ranging from 3 X 1O 1:L M0 to 1.2 X 1O 12 M0. The galaxy mass 
increases gradually from the lightest to the darkest color. The black dashed line is for the pure HBI case, where there is no stirring. The 
run is halted at early times if a cooling catastrophe occurs. 



ies and for equally spaced masses ranging from 3 x IO^Mq 
to 1.2 x 1O 12 M . The galaxy mass increases gradually from 
the lightest to the darkest color. Thick red curves are the 
corresponding lines for 200 galaxies. The black dashed line 
is for the pure HBI case. As expected, when stirring is weak, 
the HBI prevails, leading to a systematic tangential bias in 
the orientation of magnetic fields. This insulates the core 
against thermal conduction, leading to a cooling catastro- 
phe. The fields become more tangential with time and the 
cluster eventually suffers from overcooling. On the other 
hand, for increasingly vigorous stirring (i.e., increasing the 
individual masses or number of galaxies), the field becomes 
increasingly isotropic, and a cooling catastrophe is averted. 

These results are consistent with the driven turbulence 
simulations in RO1(0, and can be broadly understood in 
terms of the simple Froude/Richardson number criterion 
outlined in Ej2j The main difference in the more realistic 
scenario we present here is that the discrete nature of the 
stirrers and resonant excitation process introduces greater 
stochasticity and time-dependent fluctuations in the veloc- 
ity field and magnetic anisotropy (e.g., compare the smooth 
curves in Fig 3 & 4 of RO10 with their noisier equivalents 
Fig [3] & [4] of this paper). But the main physical conclusions 
are unchanged. It is also interesting to note that while the 
velocity dispersion is only weakly dependent on the number 
of galaxies (depending more sensitively on galaxy masses), 
the magnetic anisotropy shows somewhat greater sensitivity. 
In particular, the magnetic field anisotropy cannot simply be 
predicted from the instantaneous velocity dispersion, as in a 
naive application of a Froude/Richardson criterion. We saw 



5 Although note that all but one of the simulations in RO10 were 
adiabatic simulations; by contrast, all the simulations presented 
here simultaneously include radiative cooling. 



similar behavior in RO10, where runs with similar asymp- 
totic velocity dispersions had similar velocity anisotropies, 
but markedly different magnetic anisotropies. The advected 
magnetic field is sensitive to the integrated past displace- 
ment history of a fluid element, and not merely the instan- 
taneous velocity field. 



4.3 Evolution of Gas Temperature and Entropy 

In Figure [S] we show the evolution of temperature profiles. 
This figure is for the models where heating is more efficient. 
Specifically, it corresponds to the following pairs of param- 
eters (150,1.2), (200,0.9), (200,1.2), where the first number 
in the parenthesis is the number of galaxies and the second 
is the galaxy mass in 10 12 Mq. Progressively older profiles 
correspond to systematically brighter colors. The final time 
corresponds to 5 Gyr and the curves are plotted every 0.1 
Gyr. As is clearly seen in this figure, these models do not 
lead to the cooling catastrophe. Several features are of inter- 
est. The temperature profile does not asymptote toward an 
isothermal profile, as is genericall y the case when therma l 
conduction alone offsets cooling (iBregman fc Davidl 1 19881 : 
iGuo fc Qhllioog ; IConrov fc Ostrikeill2008l ). Despite the fact 
that we have not introduced an additional source of cen- 
tral heating such as an AGN, the cluster is able to remain 
in a thermally stable CC (i.e., with a central temperature 
which is lower than at the cooling radius) state via heat 
transport from the outer heater reservoir alone. Without 
fine-tuning, this is impossible to achieve with thermal con- 
duction alone (when the cluster either becomes isothermal or 
undergoes a cooling catastrophe). Finally, the temperature 
profile is not always monotonic, but occasionally increases 
inward — a situation which is thermodynamically impossible 
if thermal conduction alone is operating. Note that these 
fluctuations are transient; such reversals are not present in 
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Figure 4. The evolution of the anisotropy in the orientation 
of magnetic fields. Vanishing j3 corresponds to isotropic fields. 
The more negative /3 becomes, the more tangential the fields are. 
Thin blue lines are for 100 galaxies and for equally spaced masses 
ranging from 3 X 10 n M Q to 1.2 X 10 12 M Q . The galaxy mass 
increases gradually from the lightest to the darkest color. Thick 
red curves are the corresponding lines for 200 galaxies. The black 
dashed line is for the pure HBI case. Runs are terminated when 
a cooling catastrophe sets in. 

the later stages of the evolution (progressively lighter blue 
curves correspond to later times). As we shall see in i )4.5l 
all of these features hint that an additional heat transport 
process is at play: turbulent heat diffusion. 

Figure shows the temperature evolution for the pa- 
rameters where the heating is least efficient. From left to 
right shown are: (100, 0.3), (50, 0.6), (50, 0.3). Here, the 
profiles are shown more frequently then in Figure 3 to bet- 
ter capture the evolution of the system just before the im- 
minent cooling catastrophe. For the same reason, we also 
extend the radial scale to smaller distances from the cluster 
center to show how the system becomes thermally unsta- 
ble. It is evident that in all three cases, the cluster quickly 
evolves toward a cooling catastrophe. In the final stages of 
the process, the cooling is so fast at the very center that the 
gas accretion accelerates so much that adiabatic compres- 
sion in the shells surrounding the center can heat the gas up 
(e.g., see last profile in the right panel). 

Finally, in Fig [7] we show the entropy profiles (where 
entropy is defined as K = feT/i 2 ' 3 ) for the strong heat- 
ing models. The central entropy grows somewhat, consistent 
with the rise in temperature, and as might be expected if 
heating by conduction and/or turbulent heat diffusion were 
taking place. However, these profiles show that turbulent 
convection/stirring is still a relatively gentle process; we do 
not see the flat isentropic central profile which might be ex- 
pected if turbulent convection were extremely efficient. In- 
stead, the fluid always remains stably stratified by entropy, 
which steadily increases outward at all times. 

As discussed above the cluster will develop a cooling 
catastrophe if the number of galaxies and/or their masses 
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Figure 8. Time until cooling catastrophe as a function of the 
number of galaxies and galaxy mass (in units of 10 12 M Q ). Con- 
tours are in Gyr. Models that correspond to 6 Gyr (the maximum 
duration of the simulations) are thermally stable. See text for de- 
tails. 

are too small. The time it takes for the cluster to reach 
this point (essentially the effective cooling time) is plotted 
in Figure H as a function of galaxy number and mass. 
The contours are plotted every Gyr. The models that 
exhibit the effective cooling time of 6 Gyr (the maximum 
simulation run time) are thermally stable. We point out 
that in practice the models that possess cooling times > 3 
to 4 Gyr could be considered stable as they are likely to 
experience cluster mergers that may reset the conditions 
in the ICM and further slow down or essentially delay the 
cooling process. In any case, as can be seen in Figure 5, a 
substantial fraction of the models shows appreciably long 
effective cooling times. As a technical note, we add that the 
reason for the lack of monotonicity in some of the contour 
lines as a function of galaxy number is that a single random 
seed was used to generate the conditions for a given number 
of galaxies and varying galaxy masses. 



4.4 Generation of vorticity and magnetic fields 

As we have previously seen, (/-modes must be excited for stir- 
ring by galaxies to excite volume-filling turbulence. These g- 
modes also induce vorticity (equation [6]) . Vorticity is there- 
fore an excellent tracer of the growth of <?-modes. We com- 
pute the evolution of vorticity in the central 100 kpc to 
assess if g— modes are indeed generated and trapped. Figure 
[9] (left panel) shows the evolution of the square of the scaled 
vorticity for the same set of parameters as in Figure that 
shows velocity dispersion and median velocity. The scaled 
vorticity is defined as tt = (A re f /v ie { )V x v, where A rc f = 50 
kpc and v rc { = 100 km s _1 are the reference lengthscale and 
velocity, respectively. Thin blue lines are for 100 galaxies and 
red ones for 200 galaxies. Galaxy gasses range from 3 x 10 11 
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Figure 5. The evolution of temperature profiles for the strong heating models. The panels correspond to the the following pairs of 
parameters: (150,1-2), (200,0.9), (200,1.2), where the first number in the parenthesis is the number of galaxies and the second is the 
galaxy mass in 10 12 Mq, from left to right respectively. The final time corresponds to 5 Gyr and the curves are plotted every 0.3 Gyr. 
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Figure 6. The evolution of temperature profiles for the weak heating models. From left to right are the results for the following sets 
of parameters: (100, 0.3), (50, 0.6), (50, 0.3), where the first number in the parenthesis is the number of galaxies and the second is the 
galaxy mass in 10 12 Mq. The curves are shown every 0.1 Gyr. 
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Figure 7. The evolution of entropy profiles for the strong heating models. From left to right are the results for the following sets of 
parameters: (150, 1.2), (200, 0.9), (200, 1.2), where the first number in the parenthesis is the number of galaxies and the second is the 
galaxy mass in 10 12 Mq. The curves are shown every 0.1 Gyr. 
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Figure 9. The evolution of normalized vorticity (left panel) and normalized magnetic pressure. See text for definition of normalization. 
The curves correspond to the same dataset as that shown in Figure [3] and the meaning of lines is the same as in that figure. 



Mq to 1.2 x 10 12 M Q and are uniformly sampled (lighter 
colors are for lighter galaxies). The black dashed line corre- 
sponds to the pure HBI case. A clear trend for the vorticity 
to increase with time is seen in this figure, suggesting that 
g— modes are present and at least partially trapped, leading 
to the volume-filling turbulence seen. 

As discussed in ^3 a growth in vorticity might also 
lead to growth in the magnetic field; the possibility that 
magnetic fields could be turbule ntly amplified in c luster s 
has been repeatedl y raise d (e.g.. iRuzmaikin et alj Jl989l); 
ISubramanian et all (|2006l ); iRvu et all (|2008T l; ICho et al.1 
(2009)). Whilst a detailed study is beyond the scope of this 
paper, we check whether these theoretical expectations are 
satisfied in our simulations. Fig. [5] shows that the magnetic 
energy density indeed grows in tandem with vorticity, with 
more vigorous stirring corresponding to greater field ampli- 
fication. However, the characteristic growth time appears to 
be somewhat longer. Note that the simulations were initial- 
ized with extremely small magnetic fields: the initial plasma 
beta Pi 3> /Sobs, where j3 b B ~ 100 is typically measured in 
the ICM. These small initial fields were for computational 
convenience (since the MHD approximation is satisfied with 
a trivially small magnetic field), and to ensure that mag- 
netic fields never become dynamically importanlQ. Hence, 
despite growing by a factor of 50, the magnetic energy den- 
sity has not yet reached its saturated state, and is not yet 
in equipartition with turbulence. Nonetheless, the turbulent 
amplification of the B-fields, which mirrors the growth of 
vorticity, is a robust result. 



4.5 Relative contribution to gas heating 

In H4.3I we noted a number of interesting features in the 
temperature profiles of our stable clusters. They remained 
stable CC clusters, neither becoming isothermal nor devel- 
oping cooling catastrophes, as clusters stabilized solely by 
thermal conduction generally do. Furthermore, the central 
temperature showed time-dependent oscillations, sometimes 
becoming hotter than gas further out. A temperature inver- 
sion would not happen if only thermal conduction was at 
play. This behooves us to take a closer look at what actually 
stabilizes the customary thermal runaway. We have already 
discussed the effect that turbulence can have on thermal 
conduction, by tangling field lines and countering the HBI. 
However, turbulence itself can be a source of heating, ei- 
ther via viscous dissipation of turbulence, or turbulence dif- 
fusion of high entropy gas into low entropy regions (e.g., 
iDennis &: Chandranl |2005l ). and references therein). Let us 
examine these in turn. 

As long as there is sufficient separation of scales that 
an inertial range can develop (such that the energy per 
unit mass per unit time e ~ v s /l is independent of scale), 
the heating rate from dissipation of turbulent motions is 
independent of the nature of viscosity. In particular, it 
is unimportant if our numerical viscosity is different from 
the actual physical viscosity in the ICM. The heating rate 
per unit volume due to dissipation of such motions is 
jDennis fc Chandranl 120051 ): 

„ Cdiss/SWt Cdi ss £/t , , 

1 - j - —. \^>) 

l Icdd 



° Thus, as least in these simulations, the HBI is stabilized by the 
stirring motions and not by magnetic tension. 



where I is the dominant lengthscale, Ut is the energy 
density in turbulence, and t c dd = l/vt is the eddy- 
turnover time on the dominant lengthscale. To estimate 
^cdd, we can note that vorticity f2 = V x vt has units 
of t~ dd , an d that our scaled vorticity in Fig [9] is f2 2 ~ 
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Figure 10. Heating to cooling ratios for the case of steady volume-filling turbulence driven by a source function (see text and RO10 
for more details). Top row: Conductive (solid line) and convective heating to cooling ratios as a function of time for the ICM within 100 
kpc (left panel) and 50 kpc (right panel) from the cluster center, respectively. Bottom row: Conductive (left panel) and convective (right 
panel) heating to cooling ratios as a function of radius. Progressively lighter blue color denotes later times. The curves are plotted every 
~ 100 Myr. 



0.1 (A re f /50 kpc) (wref /100 km s~ 



This implies 



>..!.!• - I-"' HI" ( ^| 



(19) 



Consistently, note that the vorticity in Fig [5] indeed takes 
~ 1 Gyr to rise to its asymptotic value. This implies that 
the heating time for turbulent dissipation of motions is: 



, Uthermal 

Cheat — ~T — Cdii 



M 2 toddy ~ 10 n yr (20) 

where Ut is the thermal energy density, and we have de- 
fined the turbulent Mach number M = v t /c s (note that 
our quoted velocities vt are in 3D). While there are fac- 
tors of order unity uncertainty, it is clear that the mild sub- 



sonic motions we explore are a negligible source of heating 
via viscous dissipation (and consistent with other estimates; 
iDennis fc Chandranl l|2005l )). This also implie s that dynam- 
ical friction heating due to galaxy motions (lEl-Zant et al.l 
2004 iKim et all 120051: iKiml 120071 : IConrov fc Ostrikerll2008l : 



Birnboim fc Dekell feoiO) is not the source of heating which 



averts the cooling catastrophe in these simulations. 

On the other hand, turbulent heat diffusion is not neg- 
ligible. One can estimate its contribution from a simple mix- 
ing length prescription as in equation [11] this shows that it 
can be at least comparable to and may exceed the thermal 
conduction contribution. However, the coefficient of turbu- 
lent diflusivity, ftturb ~ ul, is only approximate and subject 
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Figure 11. Same as for Fig llOl for the model of turbulence stirred by galaxy motions. Top row: Conductive (solid line) and convective 
heating (dashed line) to cooling ratios as a function of time for the ICM within 100 kpc (left panel) and 50 kpc (right panel) from the 
cluster center, respectively. Bottom row: Conductive (left panel) and convective (right panel) heating to cooling ratios as a function of 
radius. The curves are plotted every ~ 100 Myr; progressively lighter blue color denotes later times. 



to order unity corrections. Since we have full knowledge of 
the density, velocity, temperature and magnetic fields in our 
simulations, we can attempt to directly compute the heating 
contributions from thermal conduction and turbulent heat 
convection. In particular, at a radius r we can calculate the 
inward heat flux due to conduction: 

F cond = -/testes • VT), (21) 



where eg is a unit vector pointing in the direction of the 
magnetic field and k, is the Spitzer-Braginskii conduction 
coefficient given by k — 4.6 x 10" 7 T 5/2 erg s 'cm X K , as 



well as the convective heat flux (|Parrish et al.| [2008): 

F conv = —L-k B ((v)(Sn5T) + (n)(6vST)) + ((6nST6v)) 

(22) 

where (x) is the spatial average of quantity x in the shell and 
8x is the local deviation of that quantity from its average; 
generally the second term is dominant. We can then compare 
these to the total rate of energy loss within radius r due 
to radiative cooling. We can also compute the volumetric 
heating rate due to these two processes, via H con d = V • 
F con d, H conv = V ■ F conv , although these are of course much 
noisier quantities. 

It is useful perhaps to begin by considering a case where 
the properties of turbulence are well known: the 'strong' 
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driven turbulence case of RO10, which has volume filling tur- 
bulence by construction, and rms velocities of ~ 150km s _1 . 
Conductive (solid line) and convective (dashed line) heating 
to cooling ratios as a function of time for the ICM within 
100 (50) kpc from the cluster center are shown in the upper 
left (right) panel of Fig [TUJ It is clear that conduction only 
contributes ~ 50% of the heat necessary to overcome cooling 
and convective heat flow is an important part of the energy 
budget; indeed, in the central regions turbulent advection of 
heat is the dominant heating process (note that the cluster 
is not in complete equilibrium, so the sum of the two ratios 
is not necessarily unity). The convective heat flow shows 
time-dependent fluctuations, as might be expected. In the 
bottom panels, we show the volumetric heating to cooling 
ratios as a function of radius, for conduction (left panel), 
and turbulent convection (right panel). The curves are plot- 
ted every 100 Myr; progressively lighter colors denote later 
times. The divergence of heat fluxes is a much noisier quan- 
tity, as reflected in the plots. Nonetheless, it is clear from 
the plots that convective heating dominates near the center, 
whilst conductive heating dominates further out. This rem- 
iniscent of stable hybrid AGN+ conduction heating models 
|Ruszkowski fc Begelmanl |2002h where the AGN heats the 
cluster center and conduction is important further out. 

In Fig. 1111 we show the same plots, but for the case 
where turbulence is due to stirring by galaxies. All results 
presented in this figure are for the case of 200 galaxies, each 
with 9 x 10 Mq; this is stable against a cooling catastro- 
phe. As before, conduction is only a fraction ~ 30 — 50% 
of the energy budget. However, in this case convective heat- 
ing shows dramatic oscillations as a function of time; the 
amplitude of the oscillations H conv /C ~ 5 — 10 near the 
center is much larger than in the driven turbulence case 
H COIlv /C ~ 1 — 2. The reason for this is that the domi- 
nant lengthscales of motion are comparable or larger than 
the depicted radii, as might be expected if (/-modes are ex- 
cited (since most of the energy in g-modes are in the largest 
lengthscales, comparable to the trapping radius). For in- 
stance, from H4.4I a typical lengthscale on which vorticity is 
excited is A ~ Vt^l' 1 ~ 150 kpc (v t /W0 kms _1 )(p„/0.3) _1 . 
Since fluctuations in the velocity field span larger scales 
than the ones under interest, our calculation of H conv will 
show strong time dependence (however, the calculations of 
the conductive heat flux are of course still valid. Note that 
H COIlv /C + H con d/C has to be unity on average, since the 
cluster is stabilized against a cooling catastrophe). As noted 
earlier, Poisson fluctuations in the number of galaxies in the 
core will also drive time-dependent fluctuations in the ve- 
locity field. The gas is sloshing in the potential well; we ob- 
serve this directly too in the simulations, as the gas pressure 
maximum wanders in time from the center of the potential 
well. Nevertheless, despite the breakdown of equation 1221 in 
a rigorous sense, it is clear from the amplitude of fluctua- 
tions in the bottom panels of Fig [TT] that (as in the driven 
turbulence case) conductive heating increases outwards in 
radius, while convective heating is more important near the 
center. In particular, the dominance of convective heating 
near the center, and its positive and negative fluctuations, 
allow us to understand the fluctuating temperature profiles 
seen in Fig[S] Since conduction is only a part of the energy 
budget, there is no reason for the stabilized temperature 
profile to approach isothermality. Furthermore, the reason 



why the central temperature gradient can occasionally be- 
come inverted (with the center hotter than its surroundings) 
is clear: if a high-entropy fluid element is compressed at the 
center, this will result in higher central temperatures. While 
thermal conduction seeks to make the fluid isothermal (since 
heat flows down the temperature gradient), turbulent heat 
diffusion seeks to make the fluid isentropic (since heat flows 
down the entropy gradient). In this sense, the subsonic tur- 
bulence induced by galaxies results in only mild convection, 
since as seen in Fig the gas remains convectively stable 
with entropy increasing monotonically outward. Whilst we 
have not directly calculated the diffusion of metals directly, 
this also suggests that metal mixing to larger radii will be 
somewhat enhanced (so that metals will have a broader dis- 
tribution than the galaxies), but not greatly so. Indeed, a 
mixing- length the ory calculation of met al dispersal via tur- 
bulent diffusion bv lRebusco et al.l (|2005h , who assumes levels 
of turbulence very similar to those we have simulated, shows 
excellent agreement with observations. 

Could turbulent heat diffusion alone stabilize a thermal 
runaway? We tested this hypothesis by running purely hy- 
drodynamic simulations of our most vigorous stirring case 
(200 galaxies of mass 1.2 x 10 12 Mq), without thermal con- 
duction. The results are shown in Fig. 1121 The cluster rapidly 
u ndergoes a co oling catastrophe (consistent with the results 
of lKiml (|2007l )), demonstrating that thermal conduction is 
an essential element in stabilizing the cluster. Note that this 
model already represents the most extreme choice in pa- 
rameter space of the amount of stirring possible by galaxies. 
Also in line with these results, purely hydrodynamic simula- 
tions of 'sloshing lZuHone et al.l (|2010t ) show that the cooling 
catastrophe can be delayed by not disrupted. Besides pro- 
viding the dominant source of heating in the outer regions, 
thermal conduction also reduces stabilizing buoyancy forces 
(as discussed in ij2} and thus enables more rapid, efficient 
mixing and turbulent heat diffusion. Passive scalars such as 
metals a re more efficiently ad vected in the presence of con- 
duction (|Sharma et al. 2009a), and the same is likely true of 
the advection of entropy. Thus, intriguingly, whilst neither 
process alone can stabilize the cluster, the interplay between 
turbulence and conduction does permit stability: turbulence 
enables conduction, and conduction enables turbulence. 



5 CONCLUSIONS 

Using three-dimensional MHD simulations, we have studied 
the effect of anisotropic thermal conduction and stirring mo- 
tions due to galaxies orbiting in the cluster potential on the 
effective cooling rate in cluster cool cores. Such galaxies ex- 
cite mild subsonic turbulence with vt ~ 100—200 kms" 1 . We 
find that a combination of thermal conduction and turbulent 
heat transport can stabilize the cluster, for realistic parame- 
ter choices consistent with gravitational lensing observations 
of substructure in clusters. Unlike much previous work, there 
is no subgrid physics in our simulations: we do not invoke 
sub-grid prescriptions for the topology of the magnetic field 
(which affects the effective thermal conductivity), the mag- 
nitude and volume-filling factor of turbulence, which is cal- 
culated dire ctly from the gravity /hydro solver (unlike p re- 
vious work (|Ruszkowski fc Orj[201ol : iParrish et alJl201Ch in 
which volume-filling turbulence is inserted by hand), or tur- 
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Figure 12. Evolution of temperature profiles for a purely hydro- 
dynamic run (no thermal conduction) with 200 galaxies of mass 
1.2xlO 12 M ; our strongest heating model. Curves are shown ev- 
ery 0.1 Gyr. The cluster rapidly undergoes a cooling catastrophe, 
demonstrating that thermal conduction is an essential element in 
stabilizing the cluster; purely hydrodynamic turbulent heat dif- 
fusion alone is insufficient. 

bulent heat diffusion (which is directly simulated). We have 
also simulated a cluster with significantly higher central den- 
sity than in RO10, and still found it to be thermally stable. 
Other salient points include: 

• In order for galaxies to excite volume-filling turbulence, 
rather than have turbulence confined to galactic wakes, they 
must excite <?-modes, which requires that uj st i r < wgy D , 
where ujg™ is the Brunt- Vaisala, frequency appropriate 
when thermal conduction timescales are rapid (equation 
[5J. On the other hand, overwhelming the stabilizing buoy- 
ancy forces to randomize the magnetic field requires that 
o; s ti r > ljbv- These two requirements can be simultaneously 
satisfied since uj oc l~ 2 ^ 3 for Kolomogorov turbulence; hence, 
the low-frequency, large scale modes can be trapped, while 
the high-frequency, small scale modes overcome the HBI. 

• We observed strong growth in vorticity, which is a good 
tracer of the growth of g-modes. We also observed turbulent 
amplification of B-fields in tandem with vorticity. 

• Thermal conduction provided about ~ 30 — 50% of 
the heating budget, with the rest coming from turbulent 
heat diffusion. Viscous dissipation of turbulent motions (and 
hence dynamical friction heating) is negligible. Turbulent 
heat diffusion tends to be more important in the center of 
the cluster, while conduction plays a greater role further out. 
The predominance of turbulent heat diffusion in the center — 
which is powered by motions on large scales — implies that 
it exhibits oscillations about the equilibrium temperature 
profile, and can occasionally exhibit small temperature in- 
versions as high entropy fluid elements are compressed near 
the center. However, conduction plays a crucial part of the 
story; our most extreme stirring case still suffered a cool- 
ing catastrophe if thermal conduction was omitted. Besides 



supplying heat further out in the cluster, conduction also re- 
duces stabilizing buoyancy forces and enables more efficient 
turbulent heat diffusion. It appears that turbulence enables 
conduction to operate, as well as vice-versa. The details of 
the interplay between turbulence and conduction, as well as 
the diffusion of metals in our stirring simulations, are inter- 
esting topics for future work. 

In this paper, we have focussed on a time-steady source 
of turbulence — stirring by galaxy motions — but we stress 
that other intermittent sources of turbulence such as mergers 
or AGN outbursts, can also contribute. Indeed, a sudden rise 
in heat transport processes such as conduction and turbulent 
heat diffusion due to an increase in turbulence could effect a 
CC t o NCC transition [Gup fc Oh|[2009l ; iRuszkowski fc Ohl 
|2010| ; iParrish et all I2010T I. Other processes which could 
reorient field lines in galaxy cluster include rising bub- 
bles, which cou ld amplify and st r aight e n magnetic fields 
in their wake (IRuszkowski et all 120071 ; iGuo et all 120081 ; 
IBogdanovic et al.ll2009l ). In the future, observatio ns of Fara- 
day rotation by SKA (IBogdanovic et all |2010| ) or mag- 
netic draping around ga laxies orbiting the cluster center 
l|Pfrommer fc D ursi 2010), could probe the topology of mag- 
netic field lines and test these ideas. Finally, these ideas 
about the interplay between between the thermal conduc- 
tion, the HBI and turbulence in the inner regions of the 
cluster also apply with equal force to the interplay between 
conduction, the MTI and turbulence in the outer regions of 
the cluster, which we present elsewhere (Ruszkowski et al 
2010). 
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